Density-functional theory for the spin-1 bosons in a one-dimensional harmonic trap 
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We propose the density-functional theory for one-dimensional harmonically trapped spin-1 bosons 
in the ground state with repulsive density-density interaction and anti-ferromagnetic spin-exchange 
interaction. The density distributions of spin singlet paired bosons and polarized bosons with 
different total polarization for various interaction parameters are obtained by solving the Kohn- 
Sham equations which are derived based on the local density approximation and the Bethe ansatz 
exact results for homogeneous system. Non-monotonicity of the central densities is attributed to 
the competition between the density interaction and spin-exchange. The results reveal the phase 
separation of the paired and polarized bosons, the density profiles of which respectively approach 
the Tonks-Girardeau gases of Bose-Bose pairs and scalar bosons in the case of strong interaction. 
We give the R-P phase diagram at strong interaction and find the critical polarization, which paves 
the way to direct observe the exotic singlet pairing in spinor gas experimentally. 

PACS numbers: 03.75.Mn,67.85.Fg,71.15.Mb 



I. INTRODUCTION 

Spinor Bose gases and one-dimensional (ID) system 
are both the fascinating topics in the research of cold 
atoms [il-lj]- The studies of their crossing point for ID 
spinor bosons are also attractive in theories [a-UjI and 
have been realized in recent experiments [13l4l5l |. Spinor 
gases are prepared in the optical traps where the in- 
duced electric dipole moment determines the laser-atom 
interaction and involve an ensemble of Bose atoms con- 
densed in a coherent superposition of all possible hyper- 
fine states. The early experimentally achieved spinor 
gases include ^^Na ^, \l% and «^Rb [H, [H. In the 
three-dimensional (3D) case, the spin-dependent spin- 
exchange interactions are much weaker than the spin- 
independent short-range density-density interactions, for 
example, the ratio of them are C2/C0 (Na) = 0.03 |2g | 
and C2/C0 (Rb) = —0.005 [2l| respectively. The ground- 
state wavefunction is represented by a spinor wavefunc- 
tion which minimizes the free energy [22h24| | and the 
spin-exchange interactions give rise to a rich variety of 
phenomena such as spin domains [l7'|, textures |22| . spin 
mixi ng dyn amics [25-28] . and fragmentation of conden- 
sate |29l - [3l| etc. On the other hand, ID systems can be 
realized by confining the cold atoms in strong anisotropic 
traps where the motion of atoms is effectively ID [13- 
[iSl . [33 - l36| . The interaction among the atoms can be 
tuned in the whole regime of interaction strength via the 
idea of Feshbach resonance as well as the confinement- 
induced resonance (CIR) |37H39| . These experimental 
developments have provided unprecedented opportuni- 
ties for testing the theory of ID exactly solvable many- 
body models [3, IMffl ■ 

Many theories have studied the ID spinor gases. Un- 
der the mean-field theory, Zhang and You checked the 
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validity of a Gaussian ansatz for the transverse profile in 
the weak interaction regime and a Thomas-Fermi ansatz 
(TFA) in the strong interaction regime iSJ . Hao et al. Q 
modified the Gross-Pitaevskii (GP) equations based on 
the solution of Lieb-Lininger (L-L) w] model to reveal 
that the total densities of the ID spinor bosons exhibit 
the Tonks-Girardeau (TG) [43| properties of ID scalar 
bosons when density-density interaction is strong enough. 
The detailed TG and super Tonks-Girardeau gas (STG) 
properties of ID spinor bosons have been investigated 
particularly by Deuretzbacher et al. [8| and Girardeau et 
al. [9] with the method of Bose-Fermi mapping. 

If the spin-exchange interaction can be modulated to 
the order of density-density interaction, the competition 
between these two kinds of interaction must be consid- 
ered. Cao et al. find that the ID homogeneous spinor 
bosons under co = C2 can be exactly solved with Bethe 
ansatz (B A) method . Essler et al. show its low-energy 
degrees of freedom are equivalent to a spin-charge sepa- 
ration theory of the U(l) Tomonaga-Luttinger liquid de- 
scribing the charge sector and the 0(3) nonlinear a model 
describing the spin sector ll(j| . By means of the thermo- 
dynamical Bethe ansatz (TBA) method [111 . Il2| , Lee et 
al. and Kuhn et al. give the ground state phase diagram 
and investigate the universal thermodynamics and quan- 
tum criticality of the trapped ID spinor bosons for the 
strong interaction situation. 

So far, a method is not available for the ID trapped 
spin-1 bosons in the entire region of interaction from 
weak to strong. In this paper, we develop the Hohenberg- 
Kohn-Sham density functional theory (DFT) to inves- 
tigate the ground-state properties of ID harmonically 
trapped spin-1 bosons. DFT has been widely used for 
treatin g el ectron systems with long-range Coulomb inter- 
action [4^, [43] . It also has been successfully generalized 
to cold atom systems with short-range contact interac- 
tion [50 53]. For ID cold-atom systems, the method of 
DFT based on BA results has been developed to solve 
the ground state problem of bosons |5ll . |54| . |55| , fermions 



|56l l57j and Bose-Fermi mixtures [58| . 

Here we apply this method to study how the ground 
state of ID trapped spin-1 bosons evolves along with the 
interaction strength from weak to strong. We derive the 
Kohn-Sham (KS) equations by combing the BA solutions 
and Local Density Approximation (LDA). The ground 
state densities and energies for different interactions are 
obtained by solving these equations iteratively. The pa- 
per is organized as follows. We introduce our theory in 
Sec. II and show the numerical results in Sec. III. The 
theory part includes the model, the BA equation for ho- 
mogeneous system and the KS equations for trapped sys- 
tem. In the results part, we first show the case that all 
bosons are fully paired and then for the partially paired 
case. 



II. THEORY 



P — Ni/N. The number of particle I is A'"i = S^ where 
all atoms have parallel spin forming the ferromagnetic 
phase. The number of particle II is 2N2 where N2 pairs 
of atoms are formed between two spin states s = ±1 or 
between two spin states s = 0. Here we have supposed 
that N^ ^ N^ and that Nq is even. The sign of c^^ 
determines that the interactions between the bosons are 
repulsive or attractive, while the sign of C2^ determines 
that the spin exchange interaction in the pairs are ferro- 
magnetic or anti-ferromagnetic. 



B. Bethe ansatz equations for homogeneous system 

With the particles confined in a finite ID tube with 
length L instead of a harmonic trap as in the Hamiltonian 
(H)), the system has been exactly solved by Cao et al. jTI 
with BA method under a special condition 
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We consider N spin-1 bosons of mass m confined in 
extremely anisotropic crossed optical dipole traps with 
delta-function type density-density interaction and spin 
exchange interaction between atoms. The trap is charac- 
terized by the radial and axial angular frequencies ujj_ 
and uj with corresponding harmonic oscillator lengths 
a± — y''h/muj± and a = ^filmui respectively. When lo_\_ 
^ oj, the radial Thomas- Fermi radius is small enough 
such that only axial spin domains could form [IJ, [l^ . In 
first quantized form, the Hamiltonian for ID spin-1 gas 
can be written as 



i.e., the repulsive density-density interaction equals the 
antiferromagnetic spin-exchange interaction. They gave 
the BA equations 
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where Sij are spin-1 operators, cj^ and c\^ are ID in- 
teraction parameters which can be expressed through in- 



teraction parameters gV^ in total spin S* = 0, 2 channels 
2gi^) /3 and c^^ = 
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s-wave scattering length a^ according to g^ 
2h^a%'^ /ma\ (l - Aa|^/ai) with constant A = 1.0326 
[37[. Thus cj^ and c\^ may be manipulated in wider 
range comparing with 3D spinor system. 

The number of atoms A^^, N^ and N^ corresponding 
to spin states s = +1,0,-1 are not conserved because 
the scattering between two atoms of spin s = ±1 can 
produce two atoms of spin s = and vice versa, whereas 
the total number of atoms N = N+ + A'o + A^- and total 
spin in the z component 5*^ — N+ — N- are conserved 
yet. Therefore we may consider the system is composed 
of two parts of atoms, particle I and particle II with total 
particle numbers A^ = A^i -f 2A2 and total polarization 



with e„ (x) = [x + inc') / [x — inc'), d = c/4 and c = 
2mg/h^. {ki} in the equations is the set of quasi- 
momentum and {Aq} is the set of the spin rapidity. The 
ground state energy of the system is 



E 



2m 



N 



H^^ 



(4) 



i=l 



For positive c\^ , the spin-exchange interaction is an- 
tiferromagnetic. 2A2 particles form spin singlet bound 
pairs between two spin s = ±1 atoms or between two 
spin s = atoms, whereas A^i particles are polarized. 
The equations ^ and ^ have A^i real solutions for ki 
{i — 1, • • • A^i) and A2 pairs conjugate complex solutions 
or string solutions for ka and A^ (a = 1, • • • , 2A^2), i-e. 



ka — A; ± ic , 
Aq = A; ± ic , 



(5) 



with A; (/ = 1, • • • , A^2) real numbers [3, [TT[. Inserting 
([5]) into ^ ^ and adopting the thermodynamic limit, 
i.e., the length of gas L — >■ cx) and particle numbers 

N,Ni,N2 -^ 00 with densities n = N/L, ni = Ni/L 



and 71,2 — N'2/L finite, we easily arrive at the BA integral 
equations 
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with pi (k) and p2 (A) densities of /c and A. The integral 
boundaries B and Q are determined by the conditions 
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From ^ , we get the ground state energy per unit length 
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where the binding energy of the pair in defined as 
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To solve the BA equations, we introduce the dimen- 
sionless L-L interaction parameter 7 = c/2n — m,g/h^n 
and polarization parameter p — ni/n. Let k — Bx, 
A = Qy, B = c//3i and Q = c//?2, the densities 
of quasi-momentum and spin rapidity turn out to be 
Pi (^) = 51(2^) ^nd P2 (A) = 92 (y)- The integral BA 
equations ^ and ([7]) are translated into 

2...(.)^l + ^/^.'.r(.>.(^,^) (14) 
2 fi 



^//y.2(2/)(65-60(^,| 



7^92(2/) = l + ^/_/^5i(^)(&5-fci)(^,J^) (15) 

+ ^/_/y'52(j/')(^6 + fe4-62)(|,| 



with 



(bn -bm) (.Ti,a;2) = ^n (a:^i,a;2) - bm{xi,X2) , (16) 
and 

n/4 
6„(xi,X2)= 2 , , 3- (17) 

(71/4) + [xi - X2) 

The normalization conditions (fTO| and (fTTj) are now 
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From p2|) the ground state energy per atom for the ho- 
mogeneous ID spinor gas is 

e'^°-(n,7,p) = |-^e(7,p), (20) 
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with 



and 



e (7, p) = ei (7, _p) + 62 (7, P) + 66 (7, P) (21) 
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Here (P^ and ([25)1 can be solved numerically with the 
combination of integral equations ([T4|) . ([T5|) and the nor- 
malization (ITS)) . ([T^l). The chemical potentials are taken 
as the derivatives of (l20l) as 
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where 



/i(7,p) = 3e-7|^ + (l-p)|^, (27) 



UG. UG 

/2(7,P) = 2(3e-7^-P^ 



(28) 



We see that, when p = 1, the system is in a pure 
ferromagnetic phase with spin-polarized bosons. Here 



e {l:P = 1) coincides with e (7) in the L-L model of scalar 
bosons [40[ with interaction parameter 2j. When p = 0, 
all the bosons form pairs and the system is in a pure 
antiferromagnetic phase. In the limiting case of 7 = 0, 
the system reduces to free bosons with e (7 = 0,p) — 0. 
When 7 is very strong, the integral equations (|14p and 
(ITS|l give 51 {x) w l/27r and 32 (a;) ~ I/tt, then 
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The energy per unit length in the strong interaction case 
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is composed of three parts: the energy density of A^i free 
fermions with mass m in Ref. [58|, the energy density 
of N2 free composite fermions with mass 2m, and the 
binding energy N2eb/L of the composite fermions. This 
shows that stronger interactions favors the forming of 
boson-boson pairs. The chemical potentials in this case 
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In Fig. [TJ we plot the interaction dependence of energy 
density e (7,^) for various polarization p and compare it 
with the case of scalar bosons. The upper figure is for 
p — 0, i.e., all the bosons form the pairs. Clearly we have 
ei = and 62 (7) (red dashed line), which is obtained 
numerically from the combination of ([15]), (|19p and (P^. 
increases linearly along with 7 for 7^1 and approaches 
slowly to a constant value 7r^/48 for 7 3> 1. 62 (7) and 
the inverse parabola function Cf, (7) = —7^/4 together 
dictate that e (7) is not monotonic and has a maximum 
value 0.1169 at 7 = 0.4. The lower figure is for e(7,p) 
with p = 0, 0.2, 0.4, 0.6, 0.8, 1 corresponding to the real 
lines from bottom to top. It shows that the energy func- 
tion e {"f,p) increases along with p for the same parameter 
7. For p — 1, i.e., pure ferromagnetic phase case, e (7) 
increases monotonously to constant 7r^/3 which is exactly 
the asymptotic value of scalar bosons (gray dash-dotted 
line) . 



C. Kohn-Sham equations for trapped system 

Now we consider the spin-1 bosons in a harmonic 
trap Vext (x) — mtxP'x^ 12 by means of the DFT the- 
ory based on the theorems of Hohenberg, Kohn and 
Sham. The theory enables us to deal with the energy 
and the density profile of inhomogeneous system in the 
ground state. According to the Hohenberg-Kohn the- 
orem I of DFT, the ground-state density of a bound 
system of interacting particles in some external poten- 
tial determines this potential uniquely. Denote now the 




FIG. 1: (Color Online) The function e(7,p) for the ground 
state energy of ID homogeneous spin-1 bosons. In the upper 
panel, p = 0. In the lower panel the real lines are for p = 
0, 0.2, 0.4, 0.6, 0.8, 1 from bottom to top and the dash-dotted 
line is for scalar bosons. 



space dependent densities of particle I and particle II as 
n\ (x) and 2ri2 (x), respectively. The total density is then 
n (x) = rii (x) + 2n2 (x) and the number of the particles 
are conserved separately according to 



rii (x) dx = Ni, 
112 (x) dx — N2- 



(33) 
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The ground-state energy is a functional of the densi- 
ties EQ[ni{x) Tn2{x)\. It can be decomposed as ki- 
netic energy functional of a reference noninteracting 
system T'^'^^ [ni{x) ,n2{x)], external potential energy 
functional E^xt ["-i {x) , ri2 {x)] and KS energy functional 
Eks ["-1 {x) ,n2 (x)] involving the interactions, i.e., 

Eo[ni,n2] ^ T''"^ [ni,n2]+ Eext[ni,n2] 

+EKs[ni,n2]. (35) 

We introduce two orthogonal and normalized Bose or- 
bital functionals ^i (x) and (/)2 (x) to express the densities 



ni{x) = Ni^\{x)cpi{x), (36) 

n2{x) = N2(l)*2{x)cj)2{x). (37) 

The kinetic energy functional is written as 

T'-^/[ni,n2] = -N^ j dx^l{x)^^^i{x) 

~2N2 J dxr2 (x) ^^'i'^ (^) (38) 

and the external potential energy functional is simply 

E..,[n,M=jdxV.M^)nix). (39) 

Note for each part of bosons we introduce a single or- 
bital functional assuming that the bosons are in a quasi- 
condensate state. This is different from the system of 
fermions for which the number of orbital functional is 
decided by the number of fermions. In Ref. '5S] we have 
indicated the validity of single Bose orbital functional in 
DFT. It gives the density profile of ID bosons varying 
from a standard Gaussian shape for weak interaction to 
a half-ellipse profile for strong interaction. The density 
profile for strong interaction is consistent with that of 
noninteracting fermions except the density oscillations. 
In the limit of large particle number, the difference be- 
tween the oscillating and non-oscillating profiles becomes 
imperceptible. 

Eks [ni, n2\ includes all the contribution of the inter- 
action energies. Sometimes it is partitioned as Hartree- 
Fock energy (i.e., the mean field approximation of the 
interaction energy) and exchan ge c orrelation energy [4^, 
[stI . [SSJ . Following the way in [sO, [H [H HI, we here 
treated it as an entity with the LDA, that is, the system 
can be assumed locally equilibrium at each point x in the 
external trap, with local energy per atom provided by the 
homogenous interactional system. Thus the interaction 
energy functional Eks V^i, "-2] can be formulated as 



Eks [ni,n2]~ J dxn (x) e*^°'" [ni (a;) , ^2 (x)] (40) 

where the densities ni (x) , 712 (x) are taken at point x and 
e'^"'" [ni (x) ,71-2 (x)] takes the form of Eq. ^. Note 
that both the L-L parameter 7 (x) — mg/h'^n (x) and 
the polarization parameter p (x) — ni (x) /n (x) are space 
dependent now. With the explicit form of three terms in 
([35l). the ground state energy functional are 
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(41) 



Hohenberg-Kohn theorem II guarantees that the 
ground-state density distributions are determined by 
variationally minimizing Eq with respect to ni (x) and 
ri2 (x) |48!|. That is equivalent to minimize the free- 
energy functional F — Eq — Niei — 2A^2e2 with respect 
to <pl and (l>2, where the Lagrange multipliers 61,62 are 
introduced to conserve A^i and 2A^2- Then we can get the 
KS equations 
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where the chemical potentials /ii°™, /i2°'° of the homo- 
geneous gas for densities n 1,71,2 at x are given by Eqs. 
(I25t and (HH). With the eigenvalues of (gU and ^, the 
ground state energy can be expressed as 

Eo = iViei + 2iV262 + f dxn (x) e''°'° (x) 

77l (x) III (x) dx — / 7l2 (x) /i2 (x) (ix.(44) 



With the exactly solved ghom j-^j. different 711,712 at x, 
we can solve the KS equations p2|) and (l43t together 
with the definition of orbital functional p6p and ([57)1 to 
find the density distributions 77i (x) and 7i2 (x) and then 
calculate the ground-state energy Eq from Eq. (H^ . 

We now discuss the KS equations for the limiting 
cases of weak and strong interaction. When there is 
no interaction, KS equations correctly reduce to the ID 
Schrodinger equation of simple harmonic oscillator. The 
Bose density profiles take the standard Gaussian shape 



»^i,2 (x) — — j= exp (— x^/a^) 



(45) 



When the interaction is strong, the kinetic energies in 
and ((32]) can be ignored, we have the TFA equations 
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The values of /x" and 112 are fixed by the normalization 
conditions ^^ and ((M]). Based on (gH) and (gZ]) Kuhn 
et al. give the ground state phase diagram and find that 
the singlet pairs and unpaired bosons may form a two- 
component Luttinger liquid in the strong coupling regime 
|12| . When the interaction approaches infinitely strong, 
with the limit values of chemical potential (|5n [5^ and 
the normalization conditions (1551 IM)) , we can solve (H51 





18 




16 




14 




12 


a 




^ 


10 


t+H 




o 




t/3 


8 










C 
^ 


6 



spin-1 bosons 
U=0~12.5 



fermion 




-2-10123 
X (units of a) 




-0.5 0.0 0.5 
X (units of a) 



1.5 



FIG. 2: (Color Online) Density distribution of ID trapped 
spin-1 bosons with TV = 30 and P — 0. In the upper figure, 
five black solid lines from top to bottom are respectively for 
spinor bosons with U = 0,0.1,0.5,2,12.5. Red dotted hne 
is for non-interacting Fermions with mass 2m. The lower 
figure shows the details of the peak density for even stronger 
interactions U = 12.5, 25, 150. We observe a slightly rising of 
the peak after U > 12.5. 



|47)) to obtain the following half-ellipse-like density pro- 
files 



"-1,2 [x) 



y'27Vi,2 - (a;/ai,2)' 



TTOl 9 



(48) 



where a\ = a and a2 — ^Jh/2raLu. We see that the den- 
sity of particle I is just that of A^i noninteracting har- 
monically trapped fermions with mass m. The density 
of particles II is that of N2 noninteracting fermions with 
mass 2m. It shows that for infinitely strong interaction, 
the property of particles II approaches the TG gas of 
Bose-Bose pairs. Resorting to the Bose-Fermi mapping 
method [1,|^,|43], we may map the densities of paired and 
unpaired components exactly to those of non-interacting 
Fermions 
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where Hi {x) is the Hermite polynomials. 



(49) 



III. NUMERICAL RESULTS 

We introduce the interacting parameter U — g/ahuj. 
The space-dependent Lieb-Liniger parameter is expressed 
as 7 {x) = U /an {x). For the ID spin-1 bosons with given 
N ^ P and U , we present the numerical results for densi- 
ties rii and 71.2 of particle I and II obtained by solving the 
KS equations (|42p and (|43l) with the iteration method. 
The ground state energy Eq can be obtained via the rela- 
tion (j44p . The numerical results are summarized in Figs. 
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FIG. 3: Non-monotonicity of the central density rv^' of ID 
trapped spin-1 bosons with increasing U for A'' = 30 and 
P = 0. The exists a minimum value of rv^' a,t U = 12.5 
which is shown in the inset. 



A. P = system 

We first study the system with TV = 30 and P = 0. In 
this fully paired case, the length, density and energy are 
in units of a, 1/a and hut respectively. Fig. [2]provides an 
understanding of how the density profiles change along 
with the interaction parameter for [7 = 0, 0.1, 0.5, 2, 12.5 
in the upper figure and U = 12.5, 25, 150 in the lower 
figure. With the increasing of U, the total density pro- 
file varies from a standard Gaussian-like shape charac- 
terizing the distribution of non-interacting Bose gas to 
a half-ellipse shape indicating the distribution of non- 
interacting Fermi gas. Interestingly we find that the den- 
sity profile of TV spin-1 bosons with mass m for strongly 
interacting case e.g. U — 150 overlaps that of TV/2 
noninteracting fermions with mass 2m, except the emer- 
gence of the density oscillation in the Fermionic case. It 
means that for strong interaction all atoms in the anti- 
ferromagnetic spin-1 bosons are paired with each other, 
behaving like the TG gas of Boson pairs. 
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FIG. 4: (Color Online) (a) Evolution of the ground state en- 
ergy Eo and all contributed energy terms of ID trapped spin-1 
bosons with increasing U for A'^ = 30 and P — 0. (b) the de- 
tails of energies in the weak interaction regime, (c) the details 
of external potential energy. 



of the trap. In correspondence to the non-monotonicity 
of n^'^', Eext shows non-monotonicity too, which can be 
seen clearly in Fig. SJc). Close to the critical inter- 

7mn 



in Fig. dSt), E, 



reaches its maximum E'^^^ — 59.092fia; at U — 11.5. 
But the non-monotonicity of n(°) doesn't develop visible 
effects on other energy terms. Analogous to the ground 
state energy functional e (7, p = 0) for homogeneous spin- 
1 bosons (see the upper panel of Fig. [T|), the KS energy 
Eks, i-6. the interaction energy, increases linearly in 
weak interaction regime and approaches its peak value 
^KS^ ~ 43.033?iaj at U — 1.4, corresponding to an axial 
L-L interaction parameter 7'^"^ = 0.222 (7 at the trap 
center), followed by a monotonously decreasing in strong 
interaction regime. The contributions from T^^f ^ E^xt 
and Eks together establish the ground state energy Eo 



to increase to a maximum E^ 



85.9Mhuj at U 



1.9, corresponding to 7'"' = 0.324, and then decrease 
monotonously. We notice that the value 7*^°-' = 0.324 
for Eq^'^^ are close to 7 = 0.4 for the maximal e in the 
upper panel of Fig. [TJ We understand that the compe- 
tition between the repulsion among paired bosons and 
the binding energy of the pair gives the peak of Eq at 
U = 1.9 for the same reason as in the homogeneous case, 
while that between the repulsive density-density interac- 
tion and the anti-ferromagnetic spin-exchange interaction 
gives the non-monotonicity of n^^^ at U = 12.5. 



We surprisingly notice that the density profile does 
not change with the interaction parameter monotoni- 
cally. The full tendency of the central density n'^"^ (den- 
sity at the trap center x = 0) can be seen in Fig. [3] 
and the inset shows in detail the non-monotonicity of 
n^o\ We find that central density firstly decreases to 



a minimum value n 



(0) 



4.858/a at [/ = 12.5 before 



reaching the constant central density of non-interacting 
fermions for large U. The result can be understood as 
the competition between the repulsive density-density in- 
teraction and the anti-ferromagnetic spin-exchange in- 
teraction with equal strength. The repulsive density- 
density interaction tends to increase the distance between 
bosons while the anti-ferromagnetic spin-exchange inter- 
action leads essentially attraction between s = ±1 or 
s = bosons. For U < 12.5, the prevailing repulsions 
among bosons tend to broaden them to wider space area 
and reduce the density of bosons in the trap center to a 
minimum. The anti-ferromagnetic effect is prominent for 
U > 12.5 which contracts the bosons slightly. 

Fig. |4] describes the evolution of the ground state en- 
ergy Eq and all contributed energy terms in Eq. p5p 
as a function of U. We can see the whole trend in (a) 
and the details in the mean field regime in (b). It shows 
that the kinetic energy T^'^-f decreases slowly to a con- 
stant energy as the result of interactions restraining the 
movement of atoms. However, the external potential en- 
ergy Eext increases due to the atoms occupy wider regime 



B. P / system 

For P = 1, none of the bosons can form pair and 
the system reduces to the scalar Bose gas with density- 
density interaction 2g, whose ground state properties 
have been studied by means of DFT in Ref. [53] ■ With 
the increasing of U, the density distribution and energy 
of the system approach those of a single-component TG 
gas. 

Here we study the interesting situation of a partially 
polarized spinor gas with total polarization < P < 
1. In this case the length, density and energy are in 
units of N^/^a, N^^^/a and huj respectively. In Fig. [5] 
we illustrate the density of polarized particles ni (blue 
solid lines), density of paired bosons 2n2 (red dashed 
lines) and the total density n (black dotted lines) for 
various polarization P — 0.1,0.2,0.4,0.6,0.8 (the fig- 
ures in the columns from left to right) and interaction 
U = 0.1,2,10,50 (the figures in the rows from up to 
down). Horizontal view shows that the density profile 
ni expands gradually with the increasing of P accompa- 
nied by the corresponding shrinking of 2n2- These two 
densities together result in a slightly decreasing of the 
peak density n^^' for increasing P. An analysis of ver- 
tical scope for increasing interaction parameter U (note 
that we use different vertical axis scale in the weak in- 
teraction case U = 0.1) tells us that the density 2n2 
changes smoothly from the Gaussian distribution of non- 
interacting Bose gas to the half-ellipse distribution of TG 
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FIG. 5: (Color Online) Density distributions of ID trapped spin-1 bosons at P = 0.1, 0.2, 0.4, 0.6, 0.8 (panels from left to right) 
and U = 0.1, 2, 10, 50 (panels from top to bottom). Blue solid lines denote the densities for polarized particles ni, red dashed 
lines denote the densities for paired bosons 2n2, and black dotted lines denote the total densities n. 



gas of Bose-Bose pairs. The density peak is located in 
the center of the trap for all interaction strengths. The 
behavior of density rii is, however, a little complicated. 
Though Til shows Gaussian distribution for weak inter- 
action {U = 0.1), bosons of Particle I tend to occupy 
wider space and start to be excluded from the trap cen- 
ter for intermediate interaction {U — 2,10). The single 
peak profile of rii changes into the double-peak distribu- 
tion. This reminds us the partially phase separation of 
Bose-Fermi mixture with equal mass and equal repulsive 
interaction [4J, [5^. In our case, the phase separation of 
paired and unpaired bosons occurs in the system, i.e. the 
core area of the ID spinor gas is filled with the mixture 
of paired and unpaired bosons and in the outer region we 
find cither polarized bosons for P > 0.2 or paired bosons 
for P < 0.2. For even stronger interaction {U — 50), the 
peak of rii returns back to the trap center and the den- 
sity profile ni approaches the half-ellipse distribution of 
TG gas. Nevertheless the evidence of phase separation 
becomes more prominent in the strongly interaction case. 

As a result the total density n shows the overall Gaus- 
sian distribution showing a fully mixing phase of ni and 
2712 at weak interaction {U = 0.1). When the interaction 
increases, we find a bi-modal distribution of the total 
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FIG. 6: (Color Online) Comparison of the density distribu- 
tion of ID trapped spin-1 bosons at P = 0.4 in the case of 
strong interaction. The results for both ni and 2n2 are from 
DFT for U = 50 (solid lines), TFA for U = +oo (dashed 
lines), and Ai" = 30 non-interacting fermions (dotted lines), 
respectively. Clear evidence is observed for the phase sepa- 
ration and in this case the surrounding wings are composed 
of polarized particles. The radius of vanishing density from 
DFT is obviously larger than the TFA result. 
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FIG. 7: (Color Online) Evolution of central densities n\ , Wj 
and n*"' (panels from left to right) of trapped spin-1 bosons 
with increasing interaction parameter U. Non-monotonicity 
is seen for all partially polarized cases < P < 1. 



density with particle II imposed on the top of particle I, 
i.e. the mixed core of particles I and II is surrounded by 
two wings composed of solely polarized bosons for large 
P and large U (see the several right-down panels of Fig. 
[5]). Take the P = 0.4 case as an example. We compare 
the density plots for the two components in Fig. [51 The 
solid lines are our DFT results from the iteration solution 
of the KSE Eqs. (gH) and (O, while the dashed hues are 
analytical TFA results Eq. (|48)) for infinitely strong in- 
teraction. Also shown are the densities of non-interacting 
fermions Eq. (|49p according to the Bose-Fermi mapping 
(dotted lines) for N = 30. We identify readily the density 
oscillations with 12 peaks in the polarized component rii 
in the TG limit and 9 peaks in the paired component n2 
representing the TG gas of 9 pairs of bosons. 

An interesting observation is the non-monotonicity of 
the central density value n*^°^ which already occurs in the 
fully paired case P = 0. We take a close look at the rise- 
and-fall of the peak values for polarized bosons, paired 
bosons, and the total density. Fig. [7] shows the evolution 
of central densities with increasing interaction parame- 
ter U for different P. For the central density of polarized 
particles n[ , we find it is a constant zero at P = since 
there are no Particle I in the gas and a monotonically de- 
creasing curve at P = 1 corresponding to the scalar Bose 
gas. In all partially polarized case < P < 1, the peak 
values rapidly decrease to a minimum n\ J^^^ and then 
gradually approach to a constant for strong interaction. 
The minimum is seen to move toward the stronger inter- 
action direction with increasing P and finally disappear 
for P = 1. On the contrary, the central density of paired 
particles 712 is a constant zero at P = 1 since there are 
no paired bosons in the gas and a seemingly monotoni- 
cally decreasing curve at P = corresponding to the fully 
paired bosons. Yet we know from the result of P = in 
Fig. [3] there indeed exists a minimum due to the compe- 
tition between the density interaction and spin exchange. 
The competition persists here for all partially polarized 
cases, leading to similar shallow low-lying areas, whose 
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FIG. 8: (Color Onhne) Central densities nf\ n<"' and n*") 
as a function of P for ID trapped spin-1 bosons in the strong 
interaction case. The results of DFT for U — 50 denoted 
by lines are in agreement with those in TFA for U — +00 
denoted by symbols. 



minimum is seen to move toward the weaker interaction 
direction with increasing P as shown in the middle panel 
of Fig. [71 The total central density comes from the com- 



bination of these two terms n*^°^ = rii + 2^2 • The 
right panel shows that the minimum hollow for the total 
central density remains for all partially or fully polarized 
cases < P < 1 and the tail of n*^°-' in strong interaction 
limit firstly goes up slightly followed by a nearly linear de- 
creasing for increasing P. This non-monotonicity of the 
total central density in strong interaction case {U = 50) 
is depicted in detail in Fig. [HI for increasing P, together 
with the monotonic behavior of the central density for 
polarized bosons (down- ward) and that for paired bosons 
(up- ward) . The results from DFT are compared with the 
analytical curves from the TFA. We find that in the case 
of C/ = 50 the central densities for all polarized situation 
are already very close to the limiting value in TFA, which 
are obtained from eq. pS)) as 
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We may go further into the detailed quantum phases of 
the spinor gas by defining the radii of the vanishing den- 
sities for the two components. Fig. [SI shows the phase 
diagram as a function of the global polarization P for 
U = 50, where the axial radii of the ensemble of the 
polarized and paired components are extracted from the 
numerical result of the density profiles from the KSE. 
Without loss of generality, we set the threshold values 
of the vanishing scaled density as 0.02 in the numeri- 
cal simulation. The intersection of these two radii gives 
the boundaries which divided the phase plane into three 
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FIG. 9: (Color Online) R-P phase diagram of ID trapped 
spin-1 bosons in terms of the scaled density radii K as the 
functions of total polarization P. Three quantum phases are 
identified: spin-singlet-paired phase (S), ferromagnetic spin- 
aligned phase (F) and mixed phase of the pairs and unpaired 
bosons (M). V stands for vacuum. The results of DFT for 
17 = 50 (ni, 2n2 < 0.02 are considered as vanishing scaled 
densities) and analytical results in TFA for U = -|-oo are 
compared. 



quantum phases: spin-singlet-paired bosons S, ferromag- 
netic spin-aligned bosons F, mixed phase of the pairs and 
unpaired bosons M, while V stands for the vacuum. At 
low polarization a partially polarized region forms at the 
trap center, the radius of which increases with increasing 
polarization. At a critical polarization Pc-, the partially 
polarized region extends to the edge of the cloud. When 
the polarization increases further, the edge of the cloud 
becomes fully polarized. This process can be evidently 
seen in the lowest panels in Fig. [5] Together shown in 
Fig. ini are the theoretical results in TFA. According to 
eqs. (I15|) . the radii of the vanishing densities are calcu- 
lated as 



Ri = V2NPa, 



N(l - P) 
i?2 - \l^—- -a 



(51) 



When the two radii equal to each other, we find the crit- 
ical polarization is Pc — 0.2 and the critical radius is 
Re = 0.63. The DFT results for the radii are appar- 
ently larger than the TFA estimation for both polarized 
and fully paired bosons, which can be understood easily 
from the extension of the tail of the density profile into 
outer region (see Fig. [5]). Energetically this extension 
is due to the kinetic term neglected in TFA. The criti- 
cal polarization is in agreement with the TFA crossing at 
slightly higher polarization Pc ~ 0.23 and larger radius 
Re - 0.73. 

The phenomena of phase separation at str ong interac- 
tion is consistent with the results of Ref. ^^. There 
they gave the ground state phase diagram according 
to the evolvement of Thomas-Fermi radii of ni and n2 
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FIG. 10: (Color Online) The ground state energy per atom 
of ID trapped spin-1 bosons as functions of U for different P. 
The lines are for P=0,0. 2,0.4,0. 6,0. 8,1 from bottom to top. 



along with P by the means of TBA for strong interac- 
tion case. Our results, on the other hand, are valid for 
systems in the whole interaction regime. Interestingly 
we found the double-peak structure of the density of po- 
larized bosons at intermediate interaction, which is elu- 
sive from the method in [12l |. In a seminal experiment 
on spin-imbalanced ID two-component Fermi gas tl3{ . 
phase separation of the two-spin mixture of ultracold ^Li 
atoms trapped in an array of ID tubes is reported. The 
partially polarized core is surrounded by wings which are 
composed of either a paired or a polarized Fermi gas de- 
pending on the polarization P. The pair mechanism in 
spinor gas is challenged by the density repulsive inter- 
action, which makes the phase separation more compli- 
cated. 

Finally, just as e(7,p) shows the interaction depen- 
dence for various polarizations in the homogeneous case 
and Eq in Eq. (|44l) gives the change with interaction in 
the fully paired case, we illustrate in Fig. [lU] the evo- 
lution of ground state energies Eq per atom along with 
interaction U and total polarization P in all partially 
polarized cases. The black solid line for P = here re- 
peats the result of Eq in Fig. |4l The fully polarized 
P — 1 system is equivalent to the ID trapped repul- 
sive scalar bosons in [5^ and Eq increases monotonously 
and approach the ground state energy of non-interacting 
fermion system. The energy for partially polarized sys- 
tem with < P < 1 interpolates between these two 
extremes, i.e., the positive energy terms including the 
kinetic, potential and density-density repulsive interac- 
tions together compete with the negative binding energy 
in the anti-ferromagnetically paired bosons, giving rise 
to the non-monotonically dependence of the energy on 
the interaction parameter U. In the weak interaction 
case, the repulsive interaction is prominent so that Eq in- 
creases along with U. In the strong interaction case, Eq 
decreases because the repulsive energy slowly increases 
to a constant whereas the binding energy goes downward 



11 



parabolically. Larger polarization destroys the Bose-Bose 
pairs one by one, which diminishes the effect of pairing 
binding energy and finally leads to the monotonia behav- 
ior of £;o for P =1. 



IV. CONCLUSION 

In conclusion, using DFT we study the density distri- 
bution and energy of the ID harmonically trapped spin- 
1 bosons in the ground state. We numerically solve the 
KSEs based on LDA and the solution of Bethe ansatz. 
The results show that the competition between the re- 
pulsive density-density interaction and antiferromagnetic 
spin-exchange interaction results in complicated density 
distributions and energy evolutions along with the inter- 
action parameter. We found a non-monotonic behavior 
in the central densities of both spin-singlet paired and po- 
larized bosons. Some polarized bosons are repelled out 
of the trap center in the intermediate interaction region, 
showing the double-peak structure of density profiles. 
The total density exhibits a bi-modal distribution with 
paired bosons imposed on the top of polarized bosons. 
The phenomena of phase separation occurs for strong 



interaction with the partially polarized core surrounded 
by wings which are composed of either paired bosons or 
polarized bosons depending on the polarization P. We 
give the R-P phase diagram at strong interaction and 
find that the critical polarization Pc in DFT is slightly 
larger than the TEA result. Although we treat with an 
integrable model with equal repulsive density interaction 
and antiferromagnetic spin- exchange, the results do shed 
some light on the relativistic spinor gases. We speculate 
that the new quantum phases investigated here could be 
probed in experiment by in situ imaging, analogously to 
the ID trapped Fermi gas |13| . 
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